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Abstract 

We present results from two-dimensional numerical simulations of a supersonic 
turbulent flow with parameters characteristic of the interstellar medium at the 1 kpc 
scale in the plane of the galactic disk, incorporating shear, thresholded and discrete 
star formation (SF), self-gravity, rotation and magnetic fields. A test of the model 
in the linear regime supports the results of the linear theory of Elmegreen (1991a). 
At low shear, a weak azimuthal magnetic field stabilizes the medium by opposing 
collapse of radial perturbations, while a strong field is destabilizing by preventing 
Coriolis spin- up of azimuthal perturbations (magnetic braking). At high shear, 
azimuthal perturbations are sheared into the radial direction before they have time 
to collapse, and the magnetic field becomes stabilizing again. 

In the fully nonlinear turbulent regime, while some results of the linear theory 
persist, new effects also emerge. The production of turbulent density fluctuations 
appears to be affected by the magnetic field as in the linear regime: moderate field 
strengths cause a decrease in the time-integrated star formation rate, while larger 
values cause an increase. A result not predicted by the linear theory is that for very 
large field strengths, a decrease in the integrated SFR obtains again, indicating a 
"rigidization" of the medium due to the magnetic field. Other exclusively nonlinear 
effects are: a) Even though there is no dynamo in 2D, the simulations are able to 
maintain or increase their net magnetic energy in the presence of a seed uniform 
azimuthal component, b) A well-defined power-law magnetic spectrum and an 
inverse magnetic cascade are observed in the simulations, indicating full MHD 
turbulence. Thus, magnetic field energy is generated in regions of SF and cascades 
up to the largest scales, c) The field has a slight but noticeable tendency to be 
aligned with density features. This appears to be as much a consequence of the 
gas pushing on the magnetic field as due to constraints on gas motions because 
of the presence of the magnetic field, d) A "pressure cooker" effect is observed in 
which the magnetic field prevents HII regions from expanding freely, as in the recent 
results of Slavin & Cox (1993). e) The orientation of the large-scale azimuthal field 
appears to follow that of the large-scale Galactic shear, f) A tendency to exhibit less 
filamentary structures at stronger values of the uniform component of the magnetic 
field is present in several magnetic runs. Possible mechanisms that may lead to this 
result are discussed, g) For fiducial values of the parameters, the fiow in general 
appears to be in rough equipartition between magnetic and kinetic energy. There is 



no clear domination of either the magnetic or the inertial forces, h) A median value 
of the magnetic field strength within clouds is ~ 12/iG, while for the intercloud 
medium a value of ~ is found. Maximum contrasts of up to a factor of ~ 10 
are observed. 

Subject headings: ISM: clouds - ISM: evolution - ISM: magnetic fields - ISM: 
structure - instabilities - turbulence. 
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1. INTRODUCTION 

In Paper I of this series (Vazquez-Semadeni, Passot & Pouquet 1995; see 
also Passot, Vazquez-Semadeni &; Pouquet 1994), we presented two-dimensional 
(2D) numerical simulations attempting to model the Interstellar Medium (ISM) in 
the plane of the Galaxy at the kpc scale. The model incorporates self-gravitating 
hydrodynamics, parameterized radiative cooling, diffuse heating, and a prescription 
for modeling star formation (SF). Some of the first results obtained from the 
numerical simulations were the "slaving" of the temperature and thermal pressure 
to the density in the absence of thermal instabilities due to the short thermal time 
scales, the existence of a self-sustained cycle in which the stellar energy input to 
the ISM is enough to maintain the turbulence (with an efficiency of 0.05 %), the 
formation of clouds and cloud complexes mainly through collisions of gas streams 
(turbulent ram pressure), and a large scatter in the virial ratios of cloud energies, 
although nearly-virial clouds exhibited a tendency to live longer, explaining their 
observed overabundance. 

Among the significant shortcomings of the model of Paper I, the most notorious 
one is an excessive cloud temperature ( ^ 1000 K), and consequently, too low a 
density contrast (pmax/Pmin ~ 50). Additionally, two of the most obvious omissions 
are magnetic fields and rotation of the Galactic disk. In the present paper we 
present improvements of the model along both lines: we first introduce a different 
diffuse heating mimicking the effect of shielding against background UV radiation; 
this allows the flow to reach more realistic cloud temperatures (a few hundred K) 
and density contrasts ('^ 1000). We then extend the model by incorporating both 
magnetic fields and rotation. 

The magnetic field in the ISM is ubiquitous, although its precise dynamical 
effects (e.g. Troland 1990) and topology (e.g. Trimble 1990; Heiles et al. 1993) 
are still a matter of active research and debate. Some fundamental unanswered 
questions concerning the field are: 

1) Is the field of primordial origin or is it continually generated and dissipated in 
the Galaxy? In the latter case, what are the generation/amplification mechanisms? 
(e.g. Zweibel 1987; Wielebinsky & Krause 1993; Tajima, Cable, Shibata & Kulsrud 
1992)? 

2) What is the amplitude ratio between its uniform and fluctuating components, 
and what are the typical scales of the latter? 

3) Is the field fully turbulent, or does it consist simply of a superposition of weakly 
interacting waves? 

4) What is its relative importance in the global dynamics and energetics of the ISM 
and in cloud formation and support (e.g., Falgarone & Puget 1986, Troland 1990, 
Pudritz & Gomez de Castro 1991; Elmegreen 1991a; McKee et al. 1993)? 

5) Does the field strength correlate with gas density (Mouschovias, 1976a, b; Garciar 
Barreto et al. 1987; Crutcher, Kazes & Troland 1987; Myers & Goodman 1988a)? 

6) Does the field orientation correlate with density features (cloud shapes and 
elongations) (Goodman 1991; see also Heiles et al. 1993)? 

7) Is there a tendency for equipartition between kinetic, magnetic and gravitational 
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energies in molecular clouds (MCs) as observations seem to indicate (Myers & 
Goodman 1988a,b)? 

8) Is the Alfven speed the typical velocity of propagation of disturbances in the ISM 
(Falgarone & Puget 1986; Myers & Goodman 1988a,b)? 

Concerning rotation, it has long been known that it stabilizes the gravitational 
collapse of structures in the Galactic disk (Chandrasekhar, 1961; Toomre 1964; 
Goldreich k. Lynden-Bell 1965). Two main effects are present: the conversion of 
compressive motions into shearing ones by the Coriolis force, and the "restoring 
force" against radial motions arising from the interplay between the radial 
component of gravity and the centrifugal force. Recently, a combined instability 
analysis of a self-gravitating, rotating flow with heating and cooling and magnetic 
fields has been performed by Elmegreen (1991a, 1994). In the turbulent case, the 
main mode of cloud formation identified in Paper I — that arising from random 
turbulent compressions — may also be inhibited by rotation. 

In the present paper, we present two-dimensional numerical simulations aimed 
at investigating these problems. Our simulations arc limited by two facts. First, 
since the simulations reported are 2D, they cannot produce a dynamo, which is 
intrinsically a 3D effect. However, this has the advantage that other mechanisms 
of amplification of the magnetic field can be isolated and identified. Second, the 
numerical code used is pseudospectral, with periodic boundary conditions. The 
latter imply that the centrifugal and radial gravitational forces, which depend on 
the galactocentric distance, cannot be included in the model. However, as shown in 
§2.2, all the relevant aspects of the dynamics are included in the Coriolis force and 
the large scale shear, so that the omission turns out to be inconsequential. 

The plan of the paper is as follows: §2 presents the equations of the model 
and discusses the values of the parameters; §3 describes the results of simulations 
in the linear regime, and §4 discusses the nonlinear behavior of the model, stressing 
in particular the differences between the two regimes. In §5 we discuss the general 
behavior of a simulation with fiducial values of the parameters. Finally, §6 contains 
a brief summary and a discussion of the reaches and limitations of the current 
simulations. 



As in Paper I, we use a single-fluid approach to represent the ISM with several 
source terms in order to model radiative cooling, large-scale shear, and stellar and 
diffuse heating. The computations solve the nondimensionaUzed equations 



2. THE MODEL 



2.1. Equations 
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— = V X (u X B) - usV^B, (Id) 
ot 

V'</. = p-l, (le) 

in two dimensions [on the (x,?/)-planc with d/dz = 0] using periodic boundary 
conditions. As usual, p is the density, u is the fluid velocity, e is the internal energy 
per unit mass, P is the thermal pressure, B is the magnetic field, $7 is the angular 
velocity of the rotation, and ^ is the gravitational potential. Furthermore, we will 
frequently make use of the number density n = p/rrin, where mjj is the mass of 
the hydrogen atom. We recover identically the model of Paper I by setting B = 
and = in the above equations, although some of the parameter values may have 
changed, and the diffuse heating is different (sec below, §2.2). 

We use an ideal-gas equation of state P = (7 — l)/9e, where 7 = Cp/c„ is the ratio 
of specific heats at constant pressure and volume, respectively. We take 7 = 5/3. 
The temperature is related to the internal energy by e = c^T. The variables are 
normalized to characteristic values po,Uo,Lo, and i?o, given in Table 1. We 
refer the reader to Paper I for a thorough discussion of most of the parameters and 
model terms. Section 2.2 below describes those that are new or have changed for the 
present paper. The physical dimension corresponding to the side of the integration 
box is Lq. The nondimensional parameters resulting from the normalization are: 

1/2 

the Mach number Ma = Uq/cq, where Cq = [jkTo/mH) is the adiabatic speed of 
sound at the normalizing temperature To] the Jeans number J = L^/Lj, giving the 
number of Jeans' lengths Lj = {irc^/Gpo) in the simulation. The thermal diffusivity 
is denoted kt- 

The pseudospectral scheme used here introduces no numerical viscosity, since 
in Fourier space the MHD equations become a set of coupled ordinary differential 
equations for the Fourier amplitudes. Thus, we include explicit dissipation terms 
in the equations for u and B. A hyperviscosity scheme with a operator is used, 
which confines viscous effects to the smallest resolved scales. Indeed, for a fixed 
amount of dissipation at a given small scale, the higher power of the Laplacian allows 
for a smaller effective dissipation at the large scales. This technique is widespread 
in the fluid dynamics community (see, for example, McWilliams 1984; Babiano et 
al., 1987). Note that the diffusion coefficients for the velocity and magnetic fields 
are taken to be the same and equal to i/g. 

In addition, a mass diffusion pV'^p is included in the continuity equation in 
order to smooth out the density gradients, thus allowing the simulations to reach 
higher rms Mach numbers. The effects of this term are discussed at length in Paper 
I. 

Note that the kinematic viscosity and the mass diffusion coefficients arc chosen 
so that velocity and density discontinuities are spread out a few pixels, guaranteeing 
that they can be resolved. Also, note that heat generated by the kinetic and 
magnetic dissipations is not included in the calculations because it is negligible 
compared to the stellar and diffuse heating. 

Concerning self-gravity, we note that equation (le) is a modified Poisson 
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equation appropriate for infinite media, and its application to pseudospectral 
simulations with periodic boundary conditions has been discussed in detail by 
Alecian k. Leorat (1988). Essentially, this equation represents the self-gravity of 
the density fluctuations. 

Further details on the pseudospectral integration technique we use can be found 
in Paper I. The theoretical bases on which it rests can be found in Canuto et al. 
(1988). As is standard for turbulence simulations, the initial conditions for all 
variables are Gaussian fluctuations with random phases. The typical scale size of 
the initial density, temperature and velocity fluctuations is 1/8 of the integration 
box. The typical scale size of the initial magnetic field fluctuations is specified by 
the parameter k^, discussed in §2.3.2 below. 

A few words are in order concerning the choice of units (see Table 1). The 
unit of length Lq, equal to 1 kpc in physical units, is 27r in code units. Derived 
units which involve this length also contain factors of 27r. Also, the unit of velocity 
Uo is chosen equal to the speed of sound Cq at the unit of temperature Tq. The 
fundamental density unit in the code has units of a volume density. Finally, we note 
that the time step in a typical run varies in the range .Ol-.OOl in units of 1.3 x 10^ 
years. 

2.2. Model Terms 

In this section we first give a brief summary of the forms used for the model 
terms. Only model terms that are new or have changed from Paper I are discussed 
in detail here. A detailed description of all other terms can be found in Paper I. 
Numerical values of their associated parameters are given in Table 1. 

2.2.1. Diffuse heating Fj 
As pointed out in the Introduction, the simulations presented in Paper I had one 
important shortcoming: the densest regions that can be associated with clouds still 
had temperatures much larger (^ 1, 500 K) than actual cloud temperatures (10-100 
K) together with too low a density {pm&x ~ lOpo). This problem was mainly due 
to the form of the diffuse heating used in that paper, namely a constant, uniform 
heat source. This form of the diffuse heating produced realistic temperatures for the 
Intercloud Medium (ICM) (although see the discussion in Paper I concerning the 
uncertainties surrounding the temperature of this n ^ \ cm~'^ gas), but it caused 
the medium to be too "hard" (resistive to compression), thus both reducing the 
density contrast (/Cmax/Pmin ~ 50), and keeping the cloud temperatures too high. 
Physically, if most of the diffuse background energy is in low-energy UV photons, 
then it is known that clouds can shield their interiors from this radiation if their 
column densities are high enough (Franco & Cox 1986). Conversely, if most of the 
background energy is in energetic particles, then the heating rate per unit mass 
should be rather insensitive to the cloud density. As a compromise between these 
two possibilities, we now adopt a new diffuse heating function of the form 

rd(x,t) = ro(p/pic)-" 

where piQ is a typical density of the ICM and a is a free parameter, which in 
general we take such as to give a weak dependence of Fj on p. The factor p^Q 
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in the density is introduced so that an equihbrium temperature (between diffuse 
heating and radiative cooHng) of 10'* K is obtained at the density of the ICM. This 
form of Fd is not intended to represent any realistic physical dependence, since self- 
shielding depends on column density, which is a non-local cloud property, and thus 
extremely costly to implement in our code. Instead, the adopted form provides a 
smooth density dependence. An experiment with a threshold criterion for turning 
off the diffuse heating above a critical density was presented in Paper I, but it 
was found that it introduced spurious oscillations of the density, temperature and 
pressure distributions in clouds. The fiducial values we adopt for the parameters of 
the diffuse heating are njc = Pic/^n = 0.2 cm~'^ and a = 1/2. A discussion of the 
effects of variations in a (which is reflected on an effective polytropic exponent 7eff 
[§3.1]) on the stability of the model is given in §3.2. 

2.2.2. Stellar Heating 
We model local stellar heating due to massive stars by means of a threshold 
algorithm: a heating center a few pixels across (with a Gaussian profile) is turned 
on wherever the density exceeds a threshold value pc and V • u < 0. Once a "star" is 
turned on, it stays on for a time At = 6 x 10^ yr, typical of the lifetime of OB stars. 
Note that once a star is turned on, it remains fixed with respect to the numerical 
grid. 

As described in Paper I, the star formation rate (SFR) is computed in the 
simulations as the fraction of pixels that reach pc per unit time. 

2.2.3. Cooling A 

We use the parameterization of the radiative cooling functions of Raymond & 
Cox (1976) and Dalgarno & McCray (1972), as employed by Rosen et al. (1993) and 
Rosen & Bregman (1994): 

< T < 100 X 

AiT^ WO K <T < 2000 K 

A = ( AaT^-S 2000 K <T < 8000 K 
\^T^.867 8Q00 K <T < lO^K 
. A4r-0-65 10^ K <T <4xW'^ K 

Similarly to Paper 1, the cooling and heating rates are decreased by a factor of 7 
with respect to realistic values for numerical reasons, in order to reduce the stiffness 
of the resulting system. This, however, is not expected to affect the dynamics as 
their characteristic time scales are still much shorter than dynamical time scales 
(by factors between 10 and 10^). As discussed in Paper I, the shorter thermal 
time scales have the interesting consequence that the system is capable to reach 
thermal equilibrium between cooling and diffuse heating in times much shorter 
than the characteristic time for turbulence-induced density fluctuations to occur. 
The flow thus behaves as a polytropic gas with an effective polytropic exponent 
given by the condition of thermal equilibrium (§3.1). This behavior is preserved 
as long as the thermal rates are much faster than the dynamical rates, so that 
the equilibrium temperature can be reached before the density changes appreciably. 
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as is the case of our simulations. Even in the case of strong heating due to the 
locaUzed stellar heating, thermal equilibrium between the latter and the cooling is 
achieved virtually instantaneously. In this case a strong pressure imbalance with the 
surroundings exists, due to the pointwise nature of the stellar heating, producing 
expanding motions. Thus, the dynamical effects of the stellar heating are preserved 
in spite of the reduced thermal rates we use. 

It should also be pointed out that in our simulations the gas temperature never 
exceeds 3.5 x 10^ K, since heating due to massive-star winds and supernovae is 
not considered (sec §6.2 for a discussion of limitations of the simulations). Thus, 
the strong thermal instability above 10^ K in the cooling function is never reached 
in practice, and is included only for completeness. Work including such heating 
mechanisms is in progress (Gazol et al. 1995). 

Finally, the lower cutoff of the cooling is discussed at length in Paper I, although 
note that here we have extended it to 100 K, down from 300 K in Paper I. This limit 
still does not pose significant numerical difficulties, yet it allows the simulations to 
achieve more realistic cloud temperatures. 

2.2.4. Rotation 

We assume that our integration domain is located roughly at the solar circle, 
rotating about the Galactic center at a radius Rq and with angular speed Oq, with 
the X- and y-coordinates in the model corresponding to the azimuthal and radial 
directions, respectively. In this rotating frame one must consider the Coriolis force 
—2Q,o X u, the centrifugal force Vt^Rey and the radial gravitational force gij of 
the total mass contained within radius R. In the actual Galactic disk, the latter 
two forces exactly balance each other at every radius R, and gij = —^'^{R)Rey. 
This balance is at the origin of the differential rotation of the Galactic disk and its 
associated shear. In a rotating frame characterized by a single angular velocity fio, 
the total fluid velocity is the sum of the shear Ug = UsG^ and a turbulent fluctuation 
Uturb- Then, at radius R, the difference between the centrifugal force due to the 
rotation velocity of the frame and the radial gravitational force is, to first order, 
—2Q,o{d^/dR)R^Ryey, and this excess is balanced by the Coriolis force acting on 
the shear velocity, which gives Ug = —R{dQ,/dR)y. 

The simulation of this interaction of forces and the resulting shear can be 
accomplished at two equivalent levels. At a more fundamental level, the radial 
gravitational force can be included in the equations, with an explicit R- (or, 
equivalently, a y-) dependence. At a less physical, but less numerically costly level, 
the resulting shearing profile can be simply imposed on the velocity field, without 
introducing the radial centrifugal force in the equation of motion. We have chosen 
the latter option. Note, however, that due to the periodic boundary conditions 
we use, we cannot introdiicc a monotonic function of radius. Instead, we need to 
introduce a periodic function in the integration box, the simplest of which is a 
sinusoidal profile of period one of the form 

, . . . 27ry 
Ux{y) = Aosm.—. (2) 
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This is trivially accomplished by fixing the first Fourier mode of the velocity to a 
constant. We choose the shear to have period one in order to approximate in as 
large as possible a fraction of the box the monotonic character of the actual shear 
in the Galactic disk. 

It should be kept in mind that, as given by the above expression, the shear has 
the correct sign of dil/dR (i.e., Q decreasing outwards) in only half of the integration 
box (the upper and lower quarters of the box) , and the opposite sign in the central 
half of the box. This has dynamical consequences which will be discussed in §3.2. 

Toh, Ohkitani & Yamada (1991) have developed an algorithm for simulating 
incompressible 2D flows with linear shear profiles with pseudospectral techniques 
(see also Feireisen, Reynolds & Ferziger 1981 for the compressible case). In MHD, 
a special choice of gauge for the magnetic potential can be used to accommodate 
shear profiles (Brandenburg et al. 1995). Here, for simplicity, we have adopted the 
above sinusoidal velocity profile. 

2.3. Parameters 

A detailed description of the criteria used to select most parameter values is 
given in Paper I. Table 1 gives a summary of the parameters' meanings and fiducial 
values, both in physical units and in nondimensional code units. In this section we 
only discuss the values and criteria for the new parameters related to the magnetic 
field and the Coriolis force. 

2.3.1 Rotation 

We adopt a rotational speed of 250 km s~^ (sec, e.g.. Shore 1989), implying 
O = 3.14 X 10~® yr~^, which, in units of the code becomes = 0.41. 

2.3.2. Magnetic field 
Current views picture the Galactic magnetic field along the Galactic plane in 
the solar neighborhood as having a nearly azimuthal uniform component of strength 
~ 1.5/uG, and a turbulent component Bt ^ 5fiG (e.g. Rand & Kulkarni 1989, 
Wielebinski & Krause 1993). We choose a unit of the magnetic field strength such 
that va = Uo- At the unit of density for the code (po = 1 cm~^), this corresponds 
to Bo = 5/LtG. Note that, in the nondimensional code units the Alfven speed is given 
by vl = By p. 

Finally, as mentioned in §2.1, another initial parameter of the simulations is k^, 
the characteristic wavenumber of the initial magnetic field fluctuations. Wc adopt 
a fiducial value kg = 1, although the effect of varying this parameter is described 
in §4.2. 

3. LINEAR EVOLUTION 
As will be seen in §4, the fully nonlinear behavior of the model is extremely 
complex. In order to appropriately interpret the various processes at play and to 
distinguish between linear and nonlinear effects, we first discuss the linear evolution 
of the model. The discussion in this and the following sections relies heavily on 
the results of nearly 70 runs wc performed, which we used to explore the effects of 
variations in the parameters (low resolution runs, typically 128^ grid points) and 
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to analyze the strueture and evolution of the model with the fiducial values of the 
parameters (run 28, with 512^ grid points). A summary of the parameters of the 
runs referred to in this paper is given in Table 2. 

3.1. Equilibrium state 
For the linear evolution runs, the system is started from an equilibrium state 
with uniform density Uq = 1 cm~'\ velocity field given by eq. (2), uniform 
temperature and a uniform magnetic field in the x ( "azimuthal" ) direction. Several 
cases with various values of the initial uniform magnetic field have been analyzed. 
The initial value of the temperature is that corresponding to equilibrium between 
cooling and the diffuse heating, and is thus a function of the density exponent a in 
the expression for (§2.2). Indeed, setting pA = F^ gives the following equilibrium 
values for the temperature and thermal pressure of the flow: 



eq 



Aipi+° 



i/ft 



P. 



eq 



eq _ 



l7eff 



7 



7 L Ai 



(3) 



where j^s = 1 ~ (1 + Q!)/A) and Pi and Aj are respectively the exponent and 
coefficient of the temperature in the i-th range of the cooling function. Thus, the 
initial temperature listed in Table 2 for each run, if different from 10^ K, is the 
equilibrium temperature appropriate for the value of a in that run. Note that runs 
in the fully nonlinear regime, discussed in §4, always start at T = 10^ K. For small 
perturbations about this equilibrium state, the effective speed of sound Cgfi satisfies, 



in nondimensional units, c^g 



7eff^'eq/Peq = {leff/l')Te 



eq- 



3.2. Dispersion relation and stability criteria 
The original gravitational instability analysis of Jeans (1902) has been 
continually extended by a number of authors to include additional processes, such 
as rotation (Chandrasckhar 1961; Toomre 1964; Goldreich & Lynden-Bell 1965), a 
variety of energy sources and sinks (Struck-Marcell & Scalo 1984), and all of the 
above plus magnetic fields (Elmegreen 1991a, 1994; we will refer to the latter as 
E94). In particular, the system considered by Elmegreen is nearly identical to the 
system we consider in our simulations, except that our model is two-dimensional 
(note also the reversed choice of the x- and y-axes). For the combined instability, 
E94 gives dispersion relations aX t = for radial and azimuthal perturbations (see 
also Elmegreen 1991b). They read, respectively: 



LO 



R 



2'KGak 



1 2 / 7efr 2 I 2 \ 

k (— c 



K 



(4a) 



^ Note that in Paper I the quantity c^^ used in the dispersion relation for 
gravito-acoustic waves was mistakenly written without the factor I/7. This caused 
an underestimation of the computed periods, which are then in not as good an 
agreement with the periods observed in the simulations. However, the observed and 
computed periods still agree within a factor 30%. 
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.^ = 2.G..-F^^--;^. (4.) 

where ujR { resp. uja) is the growth rate in the radial (resp. azimuthal) case, k is the 
wavenumber, a is the surface density of the disk, k = 20[1 + 1/2 (R/il) dil/dR]^^'^ 
is the epicychc frequency, and va is the Alfven velocity. Note that in E94 there is 
no factor of I/7 in the second term of the r.h.s., because in that paper c is defined 
as the isothermal sound speed, whereas in the present paper it is defined as the 
adiabatic sound speed. Also, in transcribing eqs. (4) from E94, we have omitted a 
reduction factor included there to account for the finite thickness of the disk. Note 
also that in E94 7efr is stricto sensu a different quantity than that used in the present 
paper. While ours is an equilibrium 7eff (denoted 7eg), in E94 it is a perturbation 
7eff (denoted 7gfj), which is obtained from a linear stability analysis of the internal 
energy equation, treated separately (see also Elmegreen 1991b). However, it can be 
easily shown that in the limit lo ^ 0, which corresponds to the cooling rate being 
much faster than the growth rate of the perturbation, 7^^ 'y^g. This is precisely 
the case in our simulations (cf. Paper I). Thus, it is justified to use 7eff as defined 
in §3.1 in the dispersion relation (4). Incidentally, in the opposite limit, uj 00, 
7gjj 7, where 7 is the actual heat capacity ratio of the gas. In this case, the 
growth rate of the perturbation is much faster than the cooling rate, and the flow 
responds adiabatically. 

For our 2D system without shear (k = 2i}), the dispersion relations for the 
radial and azimuthal growth rates read, in nondimensionalized form, 

Xi = J^-k\^T,^ + B',)-K\ (5a) 
7 

^A-J k — Te, ^ ^2 (56) 
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Note that, in contrast to eqs. (4), the gravitational term does not contain the factor 
k, because of the strict two-dimensionality of our model. When shear is included {k 
given by the full expression above), these growth rates cannot be defined any more 
except for t ~ 0. 

We now analyze equations (5) in the absence of shear (hence k = 2il). In our 
2D system, radial perturbations are always stable at all scales when Qr = k/J > 1 
(in the nondimensionalized units of the paper, the volume density {p) is equal to 
unity), whereas in E94 the criterion reads 



\g. 



On the other hand, in the presence of magnetic fields, azimuthal perturbations 
are always unstable (albeit slowly if is weak). The unstable wavenumbers, 
obtained by solving for in eq. (5b) and requiring A^ > 0, are such that 
k < J ^/j / {jesTeq) (the Standard Jeans wavenumber), which clearly can be smaller 
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than the smallest available wavenumber in the simulation. Thus, for practical 
purposes, our simulations with magnetic field will be unstable in the azimuthal 
direction if the unstable scale is smaller than the size of the integration box, 
i.e., > (7eff7eq/7)^/^. A similar condition for radial perturbations leads to 
> [k? + + 7eff7'eq/7)^/^. If this latter condition is satisfied, the medium 
is then unstable relative to both directions. When there is no magnetic field, the 
two dispersion relations for the radial and azimuthal directions become identical, 
and the instability criterion reduces to > (k^ + 7efrTeq/7)''■''^■ 
In the presence of shear, the above analysis is valid only for t ~ 0. Confirming 
the results in E94, our numerical simulations indicate that, for low shear rates, the 
magnetic field opposes radial perturbations which tend to compress the field lines, 
while it helps azimuthal perturbations since magnetic tension opposes Coriolis spin- 
up. For high shear, the magnetic field again stabilizes azimuthal perturbations 
because they are sheared into the radial direction before they have time to collapse 
(E94). Finally, in the non- magnetic case, there is no distinction between the radial 
and azimuthal directions, as stated before. 

Table 3 shows the results of various simulations aiming at investigating the 
combined effect of shear, magnetic fields, rotation and 7eff . All runs have J = .5. 
Shown in the Table are the relevant parameters of two sets of runs, without and 
with magnetic fields. The value of n indicated in the Table is the minimum over 
the integration box, as given by the shear profile. The seventh column gives C^, an 
indicator of the linear instability of the runs, which we define as the square of the 
ratio of the Jeans number to the appropriate rotational or thermal pressure terms 
as discussed above. The eighth column gives the actual behavior of the simulation, 
denoting by "S" and "U" stable and unstable runs, respectively. A rough estimate 
of the actual growth rate of perturbations is given by the time — in units of 10^ yr 
— taken by the simulation to reach a peak density of 5 in code units. Note that in 
Table 3 all runs have Qr, = k/ J > 1, and are thus stable to radial perturbations. 
Finally, the column labeled "collapse" gives a description of the type and/or location 
of the clouds that form. 

Note that in the present work, the initial perturbations are Gaussian with 
random phases, and therefore contain both radial and azimuthal components. A first 
observation is that many of the runs require a long time before collapsing, consistent 
with the fact that the growth of the perturbations is not exponential but oscillatory, 
due to the time-dependent nature of the linearized equations (Elmegreen 1991a). 
Indeed, a succession of epochs of growth and decay is observed in our simulations 
before gravitational collapse finally occurs in the unstable runs. 

An interesting feature of the functional form we have chosen for the shear (cf. 
eq. (2)) is that it allows representation of the effects of various amounts of shear in a 
single simulation. The profile crests have zero shear, while the nodes have maximum 
shear. Recall that, as mentioned in §2.2, the shear has the same sign as the actual 
Galactic shear {dQ,/dR < 0) in only the upper and lower quarters of the integration 
box. 

We observe that all non-magnetic runs in Table 3, except run 45, should be 
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stable according to the linear criterion based on C^. However, the linear behavior 
is seen to be modified by the presence of shear, as exemplified by runs 39 and 105, 
which differ only by the presence of shear in the former. Although both should 
be stable, since < 1, run 39 actually forms clouds in a relatively short time. 
All of these runs exhibit cloud formation in the regions where (Kl/dR < (which 
minimizes k), a reflexion of the original linear criterion. Note that the run which 
exhibits cloud formation on the shortest timescales (run 45) is the only one that 
is unstable according to the linear criterion. Moreover, this run exhibits clouds 
formation without preference for specific regions in the integration box because it 
has K = 0. 

In the presence of a magnetic field, all runs in Table 3 should be unstable 
according to the linear criterion, and indeed arc. However, the stabilizing effect 
of the magnetic field in the presence of high shear (E94) can be seen in the fact 
that the simulations with shear tend to form clouds preferentially in the regions 
of zero shear. Even though the linear instability criterion does not include the 
value of the magnetic field, nevertheless it can be observed that runs with larger 
Bo tend to form clouds more rapidly (compare runs 57 and 47; and runs 56 and 
43). This result exemplifies the inhibition of Coriolis spin- up by the magnetic field 
(magnetic braking). Also, comparing runs 42 and 43, which are identical except for 
their value of 7eff, we see that the run with zero pressure gradient collapses faster, 
again a reflection of the original linear criterion. Finally, note that the runs labeled 
"azimuthal" in the "collapse" column have no shear, leading to a collapse simply 
along field lines, and forming a radially-oriented cloud. 

These results are exemplified in figs, la and lb, in which typical unstable non- 
magnetic (fig. la) and magnetic (fig. lb) runs are respectively shown at late times in 
their evolution. In fig. la, a large elongated cloud is seen to form in the upper and 
lower quarters of the domain (recall the boundary conditions arc periodic), where 
dQ,/dR < 0. In fig. lb, clouds are seen to form at the crests of the shear profile, 
where dQ/dR = 0. 

Finally, non-magnetic runs with shear (e.g. run 41) initially exhibit an expulsion 

of material from the central regions of the box, where dfl/dR > 0, indicating that 
the classical Rayleigh instability is at work there. 

4. MAGNETIC EFFECTS IN THE FULLY TURBULENT REGIME 

4.1 Influence of star formation on magnetic field dynamics 
This section is devoted to a study of the interplay between magnetic field and 
star formation (SF) in the fully turbulent regime. We will start with the influence 
that the stellar forcing has on the magnetic field dynamics. 

Since our simulations are two-dimensional, we cannot expect any dynamo 
mechanism. Growth and maintenance of magnetic fluctuations are nevertheless 
observed in our system. In order to Tindcrstand the mechanism at play, it is 
convenient to write the equation for the evolution of the quantity V = J pA'^d'^x 
where A is the vector potential defined by the relation B = V x (^e^). After splitting 
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the magnetic field B into its constant component Bq = Boe^ and its fluctuating 
part b = V X (aez), we can write an equation for the fluctuating potential 



— = (u X B) • + ?7(-l)"+iv2^a = -u • Va + (u x Bo) • + V^^a, 



where n = 1 for standard MHD and n = 4 in our simulations. Assuming the fields 
are periodic in the domain or vanish at the boundaries, one obtains for the time 
evolution of V: 



As is well known (Moffatt, 1978) when B^, = and p is constant, the quantity 
V decreases to zero, and thus magnetic fluctuations die away, whatever forcing 
is applied on the velocity field. When density fluctuations are small enough, the 
contribution of dissipation to the evolution of V is again certainly negative. In 
the general case, one cannot assert if V will decrease monotonically, although one 
expects that in general the effect of dissipation will be a cumulative decrease in 
magnetic energy. Let us now consider the case B^ / 0. When the system is forced 
by Alfven or magnetosonic waves, the first term in the right-hand side of eq. (6) 
still vanishes and thus magnetic field fluctuations cannot be maintained. This can 
be seen if one takes for a, Uy and 6p small disturbances of the equilibrium state 
a = 0, u = 0, p = 1, proportional to the eigenvectors of the linearized equations. 
When these disturbances have an oscillatory behavior (stable case), 6p and u have 
the same parity, opposite to that of a, so that /(I + Sp)au cPx = 0. When these 
disturbances grow or decay on the other hand (unstable case), a and u have the 
same parity, opposite to that of 6p, so that the first part of the integral, / au c/^x, is 
nonzero. The first term on the r.h.s. of eq. (6) is positive and can thus balance 
or even overcome dissipation. It is precisely the situation when a single cloud 
contracts or relaxes. The fields p, a, Uy are schematically depicted in fig. 2a for 
a magnetosonic wave and in fig. 2b for a cloud contraction. This mechanism is 
probably as efficient as a dynamo to sustain magnetic fluctuations, but a dynamo 
is still needed to generate the locally constant field B^ (the constant mode B^ is 
dynamically disconnected from the other modes). 

In order to test numerically the mechanism for the generation of magnetic 
fluctuations, we have performed the series of runs denoted 12, 13, 14 in Table 2 for 
which the initial fluctuations of the magnetic field are set to 6B = 0.01 and the 
constant field Bo takes respectively the values 0, 0.03 and 0.3. Figure 3 displays the 
fluctuating magnetic energy as a function of time for the three runs. It is very clear 
that we indeed have a growth of J {SB)'^cPx. when B^ ^ 0. The values of f {SB)'^(Px. 
reached during the run are also increasing with B^ ■ The fluctuating magnetic energy 




pa^(fx = 2 1 pa{u x Bo) • e^d^x - 2r? / V"(/9a) • {Va)S:ic 
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and the density fluctuations are correlated when Bq is nonzero, whereas they are 
not when Bo = (not shown). 

A few additional points are worth noting concerning the constant field B^. 
First, its presence is a regularizing factor in inviscid incompressible MHD (Bardos, 
Sulem & Sulem, 1988). Second, it also removes neutral points (probably not all 
of them if it is too weak) and thus leads to a decrease of magnetic dissipation. 
Finally, it has also been shown that the presence of a constant field Bo hinders 
small scale turbulence (Shebalin, Matthaeus Sz Montgomery, 1981). In compressible 
MHD turbulence permeated by a uniform magnetic field, the solcnoidal cnstrophy 
is much smaller than when Bo = 0, whereas its compressible counterpart or density 
fluctuations remain unaffected (A. Broc, 1993). 

In every MHD simulation, we observe a tendency for the typical scale of the 
magnetic fluctuations to grow in time. This effect can be observed in figs. 4a and 
4b, which respectively display the magnetic spectra at times i = .65 x 10* and 
t = 2.6 X 10^ yr for run 15. Note the growth of the A; = 1 mode for the by 
component. A study in three dimensions would be required to test whether this 
mechanism persists or is a consequence of the two-dimensional assumption. But let 
us mention that in the three-dimensional incompressible case, there is an inverse 
cascade of magnetic helicity, with the magnetic energy following to a lesser extent 
(Pouquct, Frisch & Leorat, 1976; Meneguzzi, Frisch &; Pouquet, 1981; Horiuchi & 
Sato, 1986, 1988). 

Figure 5 also displays spectra of the total magnetic energy together with those 
of the solenoidal and compressive components of the velocity field. Note that 
the magnetic spectrum is almost in equipartition with the solenoidal spectrum, 
as already observed in two-dimensional compressible MHD simulations (Pouquet, 
Passot Sz Leorat, 1991; Dahlburg &; Picone, 1989; Picone and Dahlburg, 1991). 
Contrary to the case of a decaying MHD compressible turbulence, we do not observe 
a systematic growth of u • B correlations. On the other hand it is quite clear 
that for finite but moderate values of Bo, this correlation oscillates noticeably, the 
extrema increasing with Bo- When Bo is zero and the wavenumber of magnetic 
fluctuations ks is initially equal to 4, the correlation actually goes steadily to zero. 
The temporary growth of u • B correlations is probably due to a transfer of magnetic 
energy created by clump contraction, to large scale Alfven waves along Bo- These 
results are markedly different when Bo is very large because of a strong confinement 
of the clouds (see below). A phenomenon which is also common to our present and 
former MHD simulations (Pouquet et al. 1991), is that the global compressibility 
of the medium (as measured by the ratio of compressible to total kinetic energy) 
increases with the intensity of the fluctuating magnetic energy. 

It is interesting to note that the topology of the magnetic field, mainly parallel 
to the local shear, is actually enforced by the shear itself. A run with the uniform 
component of the magnetic field initially perpendicular to the direction of the shear, 
evolves so that the magnetic field becomes aligned with the shear profile. 

The field topology is almost uniform in intercloud regions and very turbulent 
in cloud complexes undergoing star formation (fig. 6). The magnetic field has 
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actually an important dynamical effect in these regions, favoring dumpiness and 
fragmentation. On the other hand, at its fiducial value, it seems rather inefficient in 
intercloud regions to influence cloud motions. We indeed observe clouds propagating 
almost in any direction with respect to the magnetic field. The turbulent character of 
the magnetic field inside clouds is also possibly at the origin of the lack of correlation 
between and p in these regions, as indicated by fig. 7b which shows a scatter 
plot of the square of the magnetic field strength as a function of density for the 
complex encircled in fig. 7a. 

Finally, the typical values of the field strength within clouds is ~ 12^G, but can 
possibly reach 30//G. This can be compared to the value in the ICM of ~ 3//G. As 
opposed to quasi-stationary models, which cannot reach such high contrasts between 
cloud and intercloud field strengths (Mc Kee et al. 1993), the present turbulent 
model is capable of producing realistic values of the magnetic field strength within 
clouds. 

4.2 Influence of the magnetic field on cloud and star formation 
The second part of this section is devoted to the influence of magnetic flelds 
on the formation of density clumps (eventually leading to star formation). As in 
the linear regime, a constant field Bo hinders cloud condensations and thus star 
formation when it is small (super-Alfvenic motions) but favors it when it is large 
(sub-Alfvenic motions). In the case of even much larger values of B^, star formation 
gets inhibited again. These results are summarized in fig. 8 which displays the total 
star formation (SFR integrated over the total lapse time of the simulations) for 
runs with different values of the constant magnetic field B^ but otherwise identical 
(runs 53, 6, 19, 54, and 61 of Table 2). If the uniform field B^ is relatively weak, 
clumps get significantly disrupted by the shear - originating either from differential 
rotation or from turbulence - before they can contract along the magnetic field; 
further contraction, that has to occur across the magnetic field, is then hindered 
by magnetic pressure. On the other hand if the uniform field is relatively strong, 
it can act against the stabilizing action of the Coriolis force through the magnetic 
braking mechanism, and clump contraction can then occur along the field before 
shear can disrupt the condensation. This trend is valid even when the magnetic 
field fluctuations are at smaller scale (see runs 11, 15 and 5 in Table 2). 

For even larger values of Bo, another effect comes into play that we call the 
"pressure cooker" effect. Shells formed by stellar heating cannot expand as much 
when the magnetic field is strong, resulting in confining clouds in complexes or very 
thick filaments. Inside these complexes, star formation induces a strongly turbulent 
state where smaller roundish clouds get formed (see fig. 9 for run 61). Figure 10 also 
shows this effect by comparing three runs with different values of the initial uniform 
magnetic field component (runs 53, 19 and 61 in Table 2). However, with such an 
inhibition of shell expansion, stellar activity is globally reduced (self-propagating SF 
is not as effective); also reduced are both the compressibility and the amount of 
u • B correlations since Alfven waves cannot get efficiently excited. Now, the cloud 
complexes have a tendency to follow the shear. This "pressure cooker" effect is 
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present even at smaller values of Bo although it is less obvious. Expansion motions 
tend to amplify magnetic fluctuations perpendicular to the velocity, and as a result 
these motions get decelerated due to the influence of the tension of the amplified 
field lines. In these regions, the field lines are perpendicular to the density gradients, 
i.e. parallel to the density features (fig. 6). 

In general, one can say that the effects of a large uniform field are somewhat 
opposite to those of a large fluctuating field. 

Planar shocks propagating perpendicular to a constant magnetic field have a 
reduced compressibility, compared to the non-magnetic case, and if the prc-shock 
state is at rest, the shock must propagate at Ug > (c^ + v\y^^. If the field gets 
amplified beyond a certain value due to nonstationary expanding motions, the shock 
will not be able to propagate and will give rise to a magnetosonic wave. The shell 
will then stop to propagate as an "entity" , and to collect matter on its path, and the 
motion of the resulting cloud will become a wave propagation. In the simulations 
we observe that shells propagate at speeds of the order of 8 km s~^ while the 
typical sound speed in the ICM is 12 km s~^ and the Alfven speed ~ 30 km 
s~^. Other effects take place when the field is curved. As mentioned above, the 
field line tension tends to oppose expanding motions, while the Lorentz force, which 
is radial and oriented inward, promotes converging motions and generates blobs. 
As a result, planar structures are less likely to subsist in the magnetic case and 
we actually observe more roundish structures when the magnetic field is stronger 
and/or at smaller scales. Also, an expanding shell will break up more easily in 
presence of a magnetic field due to the Rayleigh- Taylor (and possibly the Parker) 
instability. 

Another effect of the magnetic field is linked to the scale of its fluctuations. 
Star formation is reduced when the field fluctuations are on smaller scale (compare 
runs 53, 19 and 54 with runs 11, 15 and 5 in Table 2). This effect is the result of 
the magnetic pressure acting at a comparable scale as the one of the fluid motions 
leading to clump contraction. In a run initiated with = 4 (as opposed to the 
standard case with = 1), the inverse cascade mechanism gradually enhances the 
large-scale magnetic field, and consequently star formation increases with time in 
this run. 

Star formation is also proportional to the intensity of the magnetic fluctuations 
SB, the more so when these are on large scale, because in that case, they locally act 
as a constant magnetic field (compare runs 8 and 7 with runs 19 and 5 in Table 2). 

Before closing this section, in summary we emphasize that in the turbulent runs, 
even though we partly recover the general trends of the linear stability analysis, 
it is however clear that the formation, dynamics and morphology of clouds are 
dominated by turbulent energy injected from stars, rather than by the combined 
linear instability. 



5. GENERAL BEHAVIOR OF A FULLY TURBULENT SIMULATION 

In this section we briefly describe some noteworthy features of run 28, a 
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simulation with fiducial values of the parameters, using 512 x 512 grid points (see 
Table 2). According to our choice of parameters, one pixel corresponds to a size 
of roughly 2 pc. Although our viscosity is kept at the minimum compatible with 
the stability of our numerical scheme, it damps velocity features at scales less than 
about 10 pc. 

Figures 11a, b and c respectively show the density, pressure and temperature 
for this run at time .72 x 10^ yr. Other views of this run at different times are 
given in figs. 6 and 7a. In fig. 11a, we identify structures that can be called "giant 

complexes" , with sizes of several hundreds of pc, "complexes" , with sizes ~ 100 pc, 
and "clouds" , with sizes af a few tens of pc. However, we stress that this classification 
is rather arbitrary, since in reality there is a continuum of density structures that 
can be identified by thresholding the density at continuously varying levels. This 
fact led us in Paper I to use a continuous filling factor function. For the MHD 
simulations, this point and statistical properties of the clouds will be discussed in a 
future paper. 

One of the most striking features of the density structures in the simulation 
is their amorphous character, none of the clouds being nearly circular, nor even 
elliptical. Moreover, density structures are hierarchically nested, as observed in real 
interstellar clouds (see, e.g., Scalo 1985). This phenomenon was also observed in 
non-magnetic simulations (Vazquez-Semadcni, 1994, Paper I). 

Extremely long and thin filaments of density and pressure (with low contrasts) 
are observed. The filaments are correlated with jumps of the velocity and magnetic 
field (compare figs. 11a and b with the magnetic field components shown in figs. 12a 
and 12b). This is also apparent in both the vorticity and the current (not shown). 
Similar filamentary structures obtain in the non-MHD runs of Paper I, as well as in 
the simpler non-rotating decay runs of compressible turbulence (Passot et al., 1988). 

Star- forming ("HIT") regions have sizes of a few tens of pc, thus corresponding 
to the largest observed sizes for HII regions in real galaxies. In our simulations, 
these regions cannot be much smaller than this because of the Gaussian smoothing 
used for numerical reasons (Paper I, §2.2). "HII" regions in the simulations are most 
easily observed in the pressure and temperature images, as seen in figs, lib and 11c. 
A particularly active star- forming region is seen in fig. lib near the upper-right 
corner, and in fact this region is somewhat evocative of images of the neighborhood 
of 30 Dor in the Large Magellanic Cloud. 

Global gravitational contraction of a giant complex appears to occur in the 
upper right quadrant, as can be seen in the density field for run 28 at times 
t = A2 X lO^yr (fig. 6), t = .72 x lO^yr (fig. 11a), and t = .91 x lO^yi (fig. 
7a). However, generalized collapse is halted by SF activity, which increases the 
turbulence in the highest-density regions of the complex. Also, the kinetic and 
magnetic energies in various clouds appear to be within a factor of 5 from each 
other, suggesting rough equipartition. A detailed study of the energetics of the 
clouds and complexes in run 28 will be presented in a future paper. 
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6. CONCLUSIONS 
6.1. Summary 

In this paper we have presented an extension of the model introduced in Paper I, 
now incorporating magnetic fields, Galactic disk rotation and a density-dependent 
diffuse heating which results in a softer equation of state. The behavior of the 
model was studied both in the linear and nonlinear regimes. In the linear regime, 
the model is well described by the linear theory developed by Elmegreen (1991a, 
E94; see also Elmegreen 1991b). At low shear, a weak field stabilizes the medium 
by opposing collapse of radial perturbations, while a strong field is destabilizing by 
preventing Coriolis spin- up of azimuthal perturbations (magnetic braking). At high 
shear, azimuthal perturbations are sheared into the radial direction before they have 
time to collapse, and the magnetic field becomes stabilizing again. In the absence of 
magnetic field, the problem reduces to the classical Toomre (1964) criterion, slightly 
modified for our 2D model. 

In the nonlinear regime, a variety of interesting effects are present, many of 
them unforeseeable through linear analyses: 

An amplification mechanism for the fiuctuating magnetic field has been 
identified, that allows the maintenance of magnetic energy over the long time 
evolution of the interstellar medium. It requires the presence of a constant 
component Bo and is effective for contracting or expanding motions, such as 
the expanding shells generated by stellar heating. 

The fluctuating magnetic field is generated at the scale of a cloud and is then 

transferred to larger scales. 

The magnetic field dynamics clearly does not reduce to the propagation of 
waves. The field participates to the global magneto-hydrodynamic turbulence, 
as demonstrated by the magnetic spectra, which exhibit a developed inertial 
range. 

The importance of the magnetic field on the global dynamics is complex. 
It depends both on its topology and its intensity compared to the shearing 
motions. A strong magnetic field confines matter in big complexes where 
smaller roundish clumps arc formed, globally inhibiting star formation. The 
field is turbulent inside the clouds and straight in intercloud regions. 
No unique correlation is found between the density and the magnetic field 
intensity probably due to its turbulent character. However, its orientation at 
the edge of the clouds tends to be perpendicular to the density gradient. Also, 
the magnetic field strength within clouds is roughly a factor of 4 larger than in 
the ICM with excursions of up to factors of 10. This suggests that the field in 
clouds is amplified throTigh fiux freezing by the same collisions of gas streams 
that form the clouds (Hunter et al. 1986; Elmegreen 1993; Paper I). 
The velocity of the expanding shells rapidly becomes sub-Alfvenic. 



6.2. Discussion of limitations 
Some of the most important limitations of this work are the relatively low 
effective resolution, the two-dimensionality of the simulations, and the omission of 
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supernovae. In this section we briefly discuss their possible consequences. 

As mentioned in §5, viscosity damps velocity features at scales smaller than 3-5 
pixels (~ 6-10 pc at a resolution of 512 grid points per dimension). Additionally, the 
mass diffusion smooths out density fluctuations at comparable scales. These scales 
correspond to those of sizable molecular clouds. Thus, the simulations are incable of 
resolving the structure within the clouds of those sizes, the applicability of the model 
being restricted to structures ranging from superclouds to large individual molecular 
clouds. However, the usage of a hyperviscosity scheme with a high power of the 
Laplacian guarantees that features in this range of scales are virtually unaffected 
by viscosity, and are thus fully turbulent. Moreover, note that in this paper we 
have focused on problems that do not require a very high resolution to deal with, as 
evidenced by the fact that only one run of those reported in Table 2 has a resolution 
of 512 grid points per dimension, while all others have 128. 

The two-dimensionality of the simulations has several effects. Among them 
are the absence of vortex stretching, and the modification of the gravitational 
potential leading to a force proportional to the inverse of the distance between 
clouds. Other effects are also mentioned in Vazquez-Semadeni (1994) and Paper I. 
Moreover, a dynamo effect is obviously absent from our simulations. However, we 
have identified an alternative mechanism by which expanding bubbles can amplify 
the magnetic energy through amplification of the magnetic potential, provided a 
uniform component of the field is present. 

Another important consideration is that the two-dimensionality of the 
simulations implies the existence of an inverse (from small to large scales) energy 
cascade which is not present in three dimensions (Kraichnan 1967), and therefore a 
reduced rate of energy dissipation at small scales compared to the 3D case (see, e.g., 
Lesieur 1990). Note, however, that this does not lead, for example, to excessively 
turbulent regimes. In purely-hydrodynamic incompressible turbulence with random 
extended forcing, the inversely-cascading energy organizes itself into large-scale 
vortex pairs which exist within a normally turbulent background (McWilliams 
1984). Furthermore, in the case of our compressible simulations with a pointwise 
compressible forcing, the latter (stellar heating) produces expanding shells that 
create turbulence at large scales, although without the formation of large-scale 
vortices. This process is different from an inverse cascade and must be present 
in 3D as well as in 2D, decreasing the difference between the two cases. 

Supernovae (SNe) are omitted from the present work in part because new code 
development would be required. It is not clear whether their presence will have 
(or not) a sizable effect in particular on the density and temperature of the ICM. 
Such an omission results in the absence in our models of a hot (T 3 x 10^-10^ K, 
p ~ 10~^ cm~^) gas phase. Under the classical picture of McKee k. Ostriker (1977), 
this hot gas occupies the vast majority of the volume (filling factor ~ 1), and in this 
framework our warm ICM (T ~ 10* K, p ~ 10"-^ cm~^) ought to be be essentially 
replaced by the hot gas. On the other hand, recent numerical results of Slavin & 
Cox (1992, 1993) suggest that the filling factor of the hot gas may be as low as 
20%. In this case, our results would only be modified by the inclusion of mostly 
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isolated bubbles of hot gas with a relatively small filling factor. Which view is more 
realistic is currently a matter of debate, and depends on a variety of issues, such as 
whether a supernova remnant breaks out of its high-density surroundings before it 
has had time to cool off (e.g., Mac Low k. McCray 1988; also G. Garcia-Segura, 
private communication), and whether the magnetic field forces the surrounding 
shell to rebound (Slavin k. Cox 1992). Our finding that a "pressure cooker" effect 
reduces shell expansion for our model HII regions seems to support the Slavin k. Cox 
picture, but the degree to which the low-density gas in our simulations will change 
upon introduction of SNe cannot be assessed precisely. The simulations of Rosen 
k Bregman (1994) do include SNe and contain a sizable fraction of hot gas, but 
since their SF scheme is smooth in space and time, significant differences may arise 
upon consideration of a more realistic, spatially-discrete, thresholded SF scheme like 
the one used in the present paper. Moreover, their calculations do not include the 
magnetic field. 

In any case, in spite of our neglect of SNe, it is likely that the morphology 
and energetics of the clouds will not be dramatically affected even if the hot phase 
is pervasive, as its pressure should be comparable to that of our warm ICM and 
the dynamics induced by expanding shells is already present in our simulations and 
with comparable total kinetic energies (note that the total kinetic energy inputs 
to the ISM from OB ionizing radiation, winds, and from SNe are of comparable 
magnitude (Abbott 1982; Garcfa-Segura, Mac Low and Langer 1995)). On the other 
hand, the inclusion of forcing from SNe and OB winds may result in the medium 
being even more turbulent, and the mechanism of cloud formation by turbulent 
density fluctuations exemplified in Paper I might become even more important. 
Work considering these effects is in progress (Gazol et al. 1995), in which we will 
address questions such as the filling factors of the various gas phases in the midplane 
of the Galactic disk and the longevity of the 10^ K gas. 
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FIGURE CAPTIONS 

Fig. 1. a) Contour plot of the density field of run 39, a non-magnetic run in the 

linear regime, at t = 2.47 x 10® yr. points upwards perpendicular to the plane 
of the page. A large cloud has formed in the region where dQ/dy < (the upper 
and lower quarters of the domain; recall the boundary conditions are periodic). 
Maximum density for this plot is Pmax ~ 5. After this time, this cloud still goes 
through several oscillations before collapsing, b) Contour plot of the density field 
of run 43 at t = 2.61 x 10® yr. This is a magnetic run, also in the linear regime. 
In this case, clouds form at the crests of the shear profile given by eq. (2), where 
d^/dy = Q. Maximum density for this plot IS /^max ^ 3. 

Fig. 2. Sketches of the density p {dash-dotted line), velocity u {long-dashed line) 
and magnetic potential a {solid line) for a magneto-sonic wave {top) and for a cloud 
compression {bottom) along the y-axis, perpendicular to Bo- For waves, p and u 
have the same parity, opposite to that of a, whereas for the compression a and u 
have the same parity, opposite to that of p. 

Fig. 3. Evolution of the fluctuating component (B — Bq)^ of the magnetic field for 
three runs with different initial values of the uniform component of the magnetic 
field Bo- Solid line: run 12, Bo = 0. Dotted line: run 13, Bo = .15 pG. Dashed 
line: run 14, Bo = 1.5 p,G. Both the amplitude of the fluctuations and their time 
derivative are seen to increase with Bq. 

Fig. 4. Magnetic spectra of run 15 (initial magnetic centroid ks = 4) at two 
different times: tg = 0.65 x 10® yr {top) and tg = 2.6 x 10® yr {bottom). Solid lines 
correspond to B^ and dotted lines to By. Note that the energy in both components 
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of the magnetic field is transferred to lower modes as time increases (particularly so 

for the By component). 

Fig. 5. Velocity and magnetic spectra for run 28 at t = 2.6 x 10^ yr. Solid line: 
total magnetic spectrum. Dotted line: solenoidal velocity spectrum. Dashed line: 
compressible velocity spectrum. Note that the magnetic and the solenoidal spectra 
are within factors of a few from each other at all scales. 

Fig. 6. Density and magnetic fields for run 28 at t = 0.42 X 10^. The density 
grey scale is logarithmic, and saturates at p = 40 cm~'^. The arrow at the bottom 
right indicates a field strength of 30 /LtG. Note the strong magnetic turbulence inside 
clouds and the rather smooth character of the magnetic field in the ICM. Regions 
of alignment of the magnetic field and density features can be seen in the filament 
in the upper right corner. However there are also regions where the magnetic field 
is perpendicular to the density features such as the lower portion of the same cloud 
and also the cloud near the center of the lower left quadrant. 

Fig. 7. (a) Density field for run 28 at t = 9.1 x 10^ yr showing a circular complex of 
radius 30 pixels 60pc). The density grey scale is as in fig. 6. (b) vs. p for the 
circular region of fig. 7a. Spherically symmetric compressions should give a linear 
relationship. The scatter diagram indicates that there is no preferred compression 
geometry. 

Fig. 8. Time integral of the star formation (SF) rate for various runs as a function 

of the initial value of the uniform component of the magnetic field Bo- From smaller 
to larger values of B^, the points represent runs 53, 6, 19, 54 and 61. At small 
Bo, SF is inhibited because B^ is small enough not to counteract magnetic braking, 
but is able to prevent radial collapse of sheared condensations. Intermediate values 
of Bo counteract magnetic braking and thus promote SF. Very large values of Bq 
inhibit SF again because the magnetic field rigidifies the medium. 

Fig. 9. Contour plot of the density field for run 61 with magnetic field vectors 
superimposed. The uniform magnetic field strength Bo for this run is 10/xG. At 
this value of Bo (about six times the fiducial value), the magnetic field is almost 
unperturbed by the fluid motions but the clouds exhibit more roundish shapes. 
The magnetic field fiuctuations, even though they are small compared to ? occur 
mainly inside the clouds. 

Fig. 10. Contour plots of the density field at t = 5.2 x 10^ yr for three runs with 

progressively larger values of Bo, but otherwise identical. From left to right, run 
53 {Bo = 0), runl9 {Bo = 1.5pG), and run 61 {Bo = 10/LtG). Note the tendency 
towards more roundish structures as Bo increases. 

Fig. 11. Grey scale images of (a) density, (b) pressure and (c) temperature fields of 
run 28 at t = .72 x 10^ yr. The bright spots within the large cloud complex in the 
upper right are star-forming regions, or "HII regions" . The grey scale in a) and b) is 
logarithmic. The density grey scale saturates at /> = 40 cm~"^ and the pressure scale 
saturates at P = 5 in code units (8.3 x 10"* cm~'^ K). Note the globally low pressure 
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contrast, with typical intercloud values ranging between 3,300 and 6,700 cm~^ K, 
and typical cloud values around 10,000 cm~'^ K. Exceptions are the "HII" regions, 
where the pressure reaches 1.8 x 10^ cm~'^ K, and which are strongly saturated in 
this image. In c), the grey scale is linear and spans a temperature range of to 
20,000 K. 

Fig. 12. Grey scale images of the x- (a) and y- (b) components of the magnetic 
field for run 28 at t = .72 x 10^ yr. The grey scale ranges from —7.5 to 7.5 fiG. The 
field is most turbulent in complexes and clouds. Note the discontinuities in both 
components of the field, which are coincident with the long, thin filaments observed 
in the density and pressure fields (figs. 11a and b). 
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